/*    Copyright 2010 Tobias Marschall
 *
 *    This file is part of MoSDi.
 *
 *    MoSDi is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 3 of the License, or
 *    (at your option) any later version.
 *
 *    MoSDi is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoSDi.  If not, see <http://www.gnu.org/licenses/>.
 */

package mosdi.subcommands;

import java.util.ArrayList;
import java.util.List;

import mosdi.discovery.EvaluatedPattern;
import mosdi.discovery.MotifFinder;
import mosdi.discovery.SequenceCountThresholdSearch;
import mosdi.fa.Alphabet;
import mosdi.fa.CDFA;
import mosdi.fa.DFAFactory;
import mosdi.fa.GeneralizedString;
import mosdi.fa.IIDTextModel;
import mosdi.index.SuffixTree;
import mosdi.paa.ClumpSizeCalculator;
import mosdi.util.BitArray;
import mosdi.util.Combinatorics;
import mosdi.util.FileUtils;
import mosdi.util.Iupac;
import mosdi.util.Log;
import mosdi.util.SequenceUtils;

public class DiscoveryMultSubcommand extends Subcommand {
	
	@Override
	public String usage() {
		return
		super.usage()+" [options] <sequence(s)-file> <abelian_pattern>\n" +
		"\n" +
		"Options:\n" +
		"  -F: read patterns from file\n" +
		"  -r: simultaneously consider reverse complementary motif\n" +
		"  -t <pvalue-threshold>: search for patterns whose p-value is below threshold\n" +
		"                         (default: 1.0)\n" +
		"  -R <sequence-ratio>: only return motifs that occur in at least\n" +
		"                       floor(<sequence-ratio> * #sequences) sequences.\n" +
		"  -c <max_clump_size>: size of clump size distribution (default:8)";
	}
	
	@Override
	public String description() {
		return
		"(DEPRECATED) Discovers motifs (restricted by abelian patterns) by " +
		"optimizing the number of sequences the motifs occurs in " +
		"(and not the total number of occurrences).";
	}

	@Override
	public String name() {
		return "discovery-mult";
	}

	@Override
	public int run(String[] args) {
		parseOptions(args, 2, "Frt:PR:c:");

		// Option dependencies
		// -- none --

		// Mandatory arguments
		String sequenceFile = getStringArgument(0);
		String pattern = getStringArgument(1);

		// Options
		boolean considerReverse = getBooleanOption("r", false);
		double pvalueThreshold = getRangedDoubleOption("t", 0.0, 1.0, 1.0);
		double sequenceRatioThreshold = getRangedDoubleOption("R", 0.0, 1.0, 0.0);
		int clumpSizeDistSize = getPositiveIntOption("c", 8); 

		Alphabet iupacAlphabet = Alphabet.getIupacAlphabet();
		Alphabet dnaAlphabet = Alphabet.getDnaAlphabet();
		List<String> sequences = null;
		double[] charDist = null;

		// create list of patterns
		List<String> patternList = null;
		if (optionSet.has("F")) {
			// read patterns from file
			patternList = FileUtils.readPatternFile(pattern);
		} else {
			patternList = new ArrayList<String>(1);
			patternList.add(pattern);
		}
		
		sequences = SequenceUtils.readSequences(sequenceFile);
		charDist = new IIDTextModel(dnaAlphabet, sequences).getCharacterDistribution();
		// total sequence length
		int sequenceLength = 0;
		for (String s : sequences) sequenceLength+=s.length();
		
		SuffixTree suffixTree = SuffixTree.buildSuffixTree(dnaAlphabet, sequences, considerReverse);
		BitArray[] sequenceOccurrenceAnnotation = suffixTree.calcSequenceOccurrenceAnnotation(considerReverse?2:1);
		BitArray[] generalizedAlphabet = Iupac.toGeneralizedString("ABCDGHKMNRSTVWY").getPositions();
		MotifFinder motifFinder = new MotifFinder(suffixTree, generalizedAlphabet, considerReverse);
		
		Log.startTimer();
		long n = 0;
		for (String s : patternList) {
			Log.startTimer();
			// parse abelian pattern
			int[] abelianPattern = Iupac.parseAbelianPattern(s);
			if (abelianPattern==null) {
				Log.errorln(String.format("Invalid pattern: %s",s));
				System.exit(1);
			}
			// calculate length
			int length = 0;
			for (int i : abelianPattern) length+=i;
			// calculate expectation
			StringBuilder sb = new StringBuilder();
			for (int c=0; c<abelianPattern.length; ++c) {
				for (int i=0; i<abelianPattern[c]; ++i) sb.append(iupacAlphabet.get(c));
			}
			double singleExpectation = 0.0;
			GeneralizedString p = Iupac.toGeneralizedString(sb.toString());
			singleExpectation+=p.getProbability(charDist);
			if (considerReverse) {
				p = Iupac.toGeneralizedString(Iupac.reverseComplementary(sb.toString()));
				singleExpectation+=p.getProbability(charDist);
			}
			
			double[] dist = SequenceUtils.calculateSequenceCountDistribution(sequences, singleExpectation, s.length(),3.0,false);

			double q = 0.0;
			int threshold = 1;
			for (int i=sequences.size(); i>0; --i) {
				q+=dist[i];
				if (q>pvalueThreshold) {
					threshold=i+1;
					break;
				}
			}
			threshold = Math.max(threshold, (int)(sequenceRatioThreshold*sequences.size()));
			// Log.printf(Log.Level.STANDARD, "Sequence count distribution: %s%n", Arrays.toString(dist));
			Log.printf(Log.Level.VERBOSE, "Abelian Pattern: %s,  Threshold: %d/%d%n", s, threshold, sequences.size());
			SequenceCountThresholdSearch search = new SequenceCountThresholdSearch(threshold, sequenceOccurrenceAnnotation);
			Log.startTimer();
			motifFinder.findAbelianPatternInstances(abelianPattern, search);
			List<EvaluatedPattern> l = search.getResults();
			Log.stopTimer("Search for instances of abelian pattern in suffix tree");
			double timeInstanceSearch = Log.getLastPeriodCpu();

			// calculate number of string in expanded abelian pattern
			double logNumber = Combinatorics.logMultinomial(p.length(), abelianPattern);
			
			Log.printf(Log.Level.VERBOSE, "Identified %d/%d strings%n", l.size(),(long)Math.round(Math.exp(logNumber)));
			for (EvaluatedPattern m : l) {
				Log.printf(Log.Level.VERBOSE, "%s %d%n",iupacAlphabet.buildString(m.getPattern()), m.getScore());	
			}

			Log.startTimer();
			for (EvaluatedPattern m : l) {
				Log.startTimer();
				List<GeneralizedString> genStringList = new ArrayList<GeneralizedString>(1);
				String forwardPattern = iupacAlphabet.buildString(m.getPattern()); 
				genStringList.add(Iupac.toGeneralizedString(forwardPattern));
				if (considerReverse) {
					String reversePattern = Iupac.reverseComplementary(forwardPattern);
					if (forwardPattern.compareTo(reversePattern)>0) {
						Log.println(Log.Level.VERBOSE, "forward pattern > reverse pattern --> skipping");
						Log.stopTimer("Total time for this motif");
						continue;
					}
					genStringList.add(Iupac.toGeneralizedString(reversePattern));
				}
				// create markov additive chain
				Log.startTimer();
				CDFA cdfa = DFAFactory.build(dnaAlphabet, genStringList, 50000);
				int states = cdfa.getStateCount();
				cdfa = cdfa.minimizeHopcroft();
				int statesMinimal = cdfa.getStateCount();
				
//				// verify
//				int total = 0;
//				for (String seq : sequences) {
//					if (cdfa.countMatchesDFA(seq)>0) total+=1;
//				}
//				if (total==m.getMatches()) {
//					Log.println(Log.Level.STANDARD, "Verification ok");
//				} else {
//					throw new IllegalStateException("BUG BUG BUG");
//				}

				// calculate clump size distribution
				ClumpSizeCalculator csc = new ClumpSizeCalculator(new IIDTextModel(dnaAlphabet.size(), charDist), cdfa, forwardPattern.length());
				double[] clumpSizeDist = csc.clumpSizeDistribution(clumpSizeDistSize, 1e-30);

				Log.restartTimer("calculate clump size distribution");
				double timeClumpSizeDist = Log.getLastPeriodCpu();
				// calculate expected clump size
				double expectedClumpSize = 0.0;
				for (int i=1; i<clumpSizeDist.length; ++i) {
					expectedClumpSize+=clumpSizeDist[i]*i;
				}
				
				dist = SequenceUtils.calculateSequenceCountDistribution(sequences, singleExpectation, m.getPattern().length, expectedClumpSize, false);
				double pvalue = dist[m.getScore()];
				Log.stopTimer("p-value calculation");
				double timeConvolution = Log.getLastPeriodCpu();

				Log.stopTimer("Total time for this motif");
				double timeMotif = Log.getLastPeriodCpu();

				if (pvalue>pvalueThreshold) {
					Log.printf(Log.Level.VERBOSE, "pvalue to high (%e, threshold is: %e) --> skipping motif%n",pvalue, pvalueThreshold);
				} else {
					sb = new StringBuilder();
//					if (logarithmic) {
//						sb.append(String.format(">>p_value>> %s LOG ", LogSpace.toString(pvalue)));
//					} else {
						sb.append(String.format(">>p_value_seq>> %e LIN ", pvalue));
//					}
					sb.append(String.format(">>stats>> %s %d %d %d %d %e %e ", forwardPattern, genStringList.size(), states, statesMinimal, m.getScore(), singleExpectation, expectedClumpSize));
//					sb.append(String.format(">>poisson>> %e %e ", lambda, expectedClumpSize));
					sb.append(Log.format(">>runtimes>> %t %t %t ", timeClumpSizeDist, timeConvolution, timeMotif));
//					if (pValueTable) {
//						sb.append(">>p_value_table>> ");
//						for (double d : dist) sb.append(String.format("%e ", d));
//					}
					if (sb.length()>0) Log.println(Log.Level.STANDARD, sb.toString());
				}
			}
			Log.stopTimer("Checking found patterns");
			double timePatternCheck = Log.getLastPeriodCpu();
			Log.stopTimer("Processing abelian pattern");
			double timeTotalAbelian = Log.getLastPeriodCpu();
			n+=l.size();
			Log.printf(Log.Level.STANDARD, ">>> %s %e %d %d %d %t %t %t %e%n", s, singleExpectation, threshold,l.size(), (long)Math.round(Math.exp(logNumber)), timeInstanceSearch, timePatternCheck, timeTotalAbelian, pvalueThreshold);
		}
		Log.stopTimer("Total time");
		double timeTotal = Log.getLastPeriodCpu();
		Log.printf(Log.Level.STANDARD, ">>!>total_time>>: %t %n", timeTotal);
		Log.printf(Log.Level.VERBOSE, "Total strings found: %d%n", n);
		return 0;
	}
}
